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We consider the influence of the perturbative bulk viscosity on the evolution of the Hubble pa- 
rameter in the QCD era of the early Universe. For the geometry of the Universe we assume the 
homogeneous and isotropic Friedmann-Lemaitre-Robertson- Walker metric, while the background 
matter is assumed to be characterized by barotropic equations of state, obtained from recent lattice 
, QCD simulations, and heavy-ion collisions, respectively. Taking into account a perturbative form 

^n^j ■ for the bulk viscosity coefficient, we obtain the evolution of the Hubble parameter, and we compare 

it with its evolution for an ideal (non-viscous) cosmological matter. A numerical solution for the 
viscous QCD plasma in the framework of the causal Israel-Stewart thermodynamics is also obtained. 
Both the perturbative approach and the numerical solution qualitatively agree in reproducing the 
viscous corrections to the Hubble parameter, which in the viscous case turns out to be slightly 
5-H , different as compared to the non-viscous case. Our results are strictly limited within a very narrow 

temperature- or time-interval in the QCD era, where the quark-gluon plasma is likely dominant. 

^ _ PACS numbers: 98.80.Es, 98.80.Cq, 12.38.Mh, 25.75.-q 
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The bulk viscosity is assumed to play an essential role in various eras of the early Universe [1| . Causal bulk viscous 
C^i \ thermodynamics has been extensively used for describing the dynamics of the early Universe, and in astrophysical 
applications. Due to the complicated nature of the cosmological evolution equations in the early Universe, very few 
exact solutions are known. Isotropic homogeneous Universes filled with causal viscous fluid obeying the relation 
£ ~ p s , where £ is the bulk viscosity coefficient, and p is the energy density, have been studied for special cases of the 
^ ' exponent s, s = 1 and s = 1/2, respectively, in and 0411, respectively. Arbitrary values of the exponent s are also 
^ \ considered in Refs. @, @. It has been proposed that causal bulk viscous thermodynamics can give a model for the 
phenomenological matter creation in the early Universe @, Q . 

Recent RHIC results apparently indicate that hot and dense matter has been formed in heavy- ion collisions Q, 
which likely agrees with the new state of matter predicted in lattice QCD simulations fl2l . [F3| and clearly points out 
that the bulk viscosity £ is non- negligible in regions close to the QCD critical temperature, T c . Nevertheless, there 
has been a debate on how to modify the cosmological standard model @, and the possibility of applying the QCD 
barotropic equations of state (EoS) [§] in the QCD era of the early Universe 0, [l(| . Reasons for such a "categorical 
resistance" might be the mathematical difficulties associated with the second Abel type non-homogeneous and non- 
linear differential equations describing the dynamics of causal bulk viscous models [H, EH EH • The long history of 
cosmological studies in which background matter has been modeled as an ideal (non-viscous) fluid would represent 
another reason. Obviously, the viscous characteristics of cosmological matter likely leads to important consequences 
in cosmology and astronomy [l4j . 

It is the purpose of the present work to investigate the QCD cosmological era of the evolution of the Universe. 
We assume that the Universe, described by the standard Friedmann-Lemaitre-Robertson- Walker (FLRW) cosmic 
background geometry, is filled with a relativistic viscous QCD plasma, whose bulk viscosity is supposed to be finite. The 



I. INTRODUCTION 



equation of state of the plasma is obtained from recent heavy-ion collisions and lattice QCD simulations, respectively 
[TBI . [l6| . The QCD plasma is assumed to be formed at temperatures in the range 0.2 < T < 10 GcV, or in the 
time-period 18.35 < t < 0.0073 GeV -1 , respectively. 
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To fully describe the effects of bulk viscosity on the cosmological evolution we analyze several scenarios. Firstly, 
we assume a perturbative bulk viscosity, £ = a, where a is a constant independent of T and/or p. However, this 
assumption has an insufficient thcrmodynamical motivation. Nevertheless, in the absence of a better alternative, we 
can consider this approach in order to obtain some indications for the deviation from the zero bulk viscosity case £ = 0. 
On the other hand the assumption of the constant bulk viscosity coefficient is physically motivated, if perturbative 
methods are used to describe bulk viscous processes. Thus, the T- and/or the p-dependence of £ or a can be ignored. 
A better approach could be the use of Af = 4 SYM [l7|, which is apparently valid at very high T. 

The main concern about considering a constant £ stems from the corresponding T- and/or p-dependence of the 
relaxation time r, and particle decay width T, which has to ensure that the speed of viscous pulses does not exceed the 
speed of light. For example, a T-dependence of r of the form r = (2 — ln2)/(27rT), guarantees a positive heat capacity. 
A constant relaxation time can be interpreted as indicating that the relaxation towards the thermal equilibrium is 
very slow. On the other hand, both the 1/T and the r times scale become arbitrarily small with decreasing time t. 
At T > T c , the quark masses are much smaller than T. Then the only way to construct a characteristic time scale for 
t is 1/T. Also, powers of dimensionless coupling a s can be multiplied by 1/T. Once again, the perturbative methods 
are consistent with such an approximation. 

Secondly, we consider numerical solutions for the evolution equation in the case of finite viscous background matter 
characterized by the Israel-Stewart causal thermodynamics. Barotropic EoS for both £ and r are adopted from recent 
lattice QCD simulations and quasi-particle (QP) effective models respectively. The numerical methods are very 
much sensitive to the boundary conditions. To keep this sensitivity as small as possible, we use the results obtained 
from the perturbative treatment. 

The present paper is organized as follows. The QCD plasma model is described in Section [TTJ The perturbative 
treatment and the numerical solution are considered in Section IIIII We discuss and conclude our results on the 
evolution of the Hubble parameter in the QCD era in Section [IVJ 



II. GEOMETRY, FIELD EQUATIONS, AND CONSEQUENCES 

In spherical coordinates (t, r, 0, </>), the line element of a homogeneous and isotropic flat Universe can be represented 
in the standard FLRW form as 

ds 2 = dt 2 - a 2 (t) [dr 2 + r 2 (d6 2 + sin 2 0d<p 2 )] , (1) 

where a is the scale factor. The Hubble parameter is defined as H — a/ a. For a vanishing cosmological constant A, 
the dynamics of the Universe is described by the Einstein gravitational field equations, given by 

Rik - \gikR = 8tt G T ik . (2) 

The energy-momentum tensor of the bulk viscous cosmological fluid filling the very early Universe is given by 

Tt = (p + p + H)u lU k -(p + U)St (3) 

where i, k takes the values 0, 1, 2, 3, p is the mass density, p the thermodynamic pressure, n the bulk viscous pressure, 
and Ui is the four velocity, satisfying the condition uiu 1 = 1. The particle and entropy fluxes arc defined according to 
N l = nu l and S* = sN l - (tW 2 /2£T) u\ where n is the number density, s the specific entropy, T > the temperature, 
£ the bulk viscosity coefficient, and r > the relaxation coefficient for transient bulk viscous effect (i.e. the relaxation 
time) , respectively. In the following we shall also suppose that the energy- momentum tensor of the cosmological fluid 
is conserved, that is T£ k = 0. 

The bulk viscous effects can be generally described by means of an effective pressure n, formally included in the 
effective thermodynamic pressure p e f / = p + n . Then in the comoving frame the energy momentum tensor has 
the components Tq = p, = T 2 = T| = —p e ff- For the line element given by Eq. (fTJ), the Einstein field equations 
read 

a J S 

d 4-7T „ , 

- = — —G (3p eff + p), (5) 

where the dot denotes the derivative with respect to the time t, and G is the gravitational constant. The energy 
density of the cosmic matter fulfills the conservation law: 



p + 3H ( Peff + p) = 0. 



(6) 
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In presence of bulk viscous stress II, the effective thermodynamic pressure term becomes p e ff = p + H. Then Eq. ([6]) 
can be written as 

p + 3H (p + p) = -3ELH. (7) 

For the evolution of the bulk viscous pressure we adopt the causal evolution equation , obtained in the simplest 
way (linear in II) to satisfy the iJ-theorem (i.e., for the entropy production to be non- negative, = II 2 /£T > 
|2ll. |22J| ) . According to the causal relativistic Israel- Stewart theory, the evolution equation of the bulk viscous pressure 
reads JTlJ 

ril + n = -3£ff - \tU (zH + I - | - ^ . (8) 

In order to close the system of equations (HJ), ([7]) and flSJ) we have to add the equations of state for p and T. The equation 
of state, the temperature and the bulk viscosity of the quark-gluon plasma (QGP), can be determined approximately 
at high temperatures from recent lattice QCD calculations as [HI, [l6| 

P = cjp, T = (3p r , e = "p+— T c 4 , (9) 



1 97 2 - 24 7 + 16 
9loq 7 — 1 



(10) 



where w = (7 — 1), 7 ~ 1.183, r ~ 0.213, /3 ~ 0.718, and uio ~ 0.5 — 1.5 GeV, respectively. In the following we assume 
that ap >> 9/ojqT^, and therefore we take £ ~ ap. In order to close the system of the cosmological equations, we 
have also to give the expression of the relaxation time r, for which we adopt the expression [Tl| . 

T = tp- 1 ~a. (11) 

Eqs. ^ are standard in the study of the viscous cosmological models, whereas the equation for r is a simple 
procedure to ensure that the speed of viscous pulses does not exceed the speed of light. Eq. (fTTj) implies that the 
relaxation time in our treatment is constant, but strongly depends on the EoS. These equations are without sufficient 
thermodynamical motivation, but in the absence of better alternatives, we shall follow the practice of adopting them 
in the hope that they will at least provide some indication of the range of bulk viscous effects. The temperature law 
is the simplest law guaranteeing positive heat capacity. 



III. PERTURB ATIVE AND NUMERICAL SOLUTIONS OF THE FIELD EQUATIONS 

For a viscous FLRW Universe filled with a viscous quark-gluon plasma (QGP), the evolution equation of the Hubble 
parameter is given by 

aHH + -[1 + (1 - r)-f]aH 2 H + HH - (1 + r)aH 2 + ^(7 - 2)aH i + ^H 3 = 0. (12) 

The parameters, a, r and 7 are the coefficients and the exponents of the barotropic QGP EoS, discussed later in Eqs. 
([9j)-(fT0"|) 0- In the limit of a vanishing a, a non-viscous evolution equation can be obtained from Eq. (fT2|) . Then, 

H (t)Ho(t) + ^H 3 (t) = 0, (13) 

where Ho(t) is the Hubble parameter corresponding to the non-viscous background matter. Eq. (|13p can easily be 
solved by 

H a (t) = ^-t-\ (14) 
37 

The mathematical solution Hq = of Eq. (|13|) . which can be interpreted physically as describing a static Universe 
with a vacuum background geometry, has been ignored, as irrelevant to the present case. The second solution refers 
to a dynamic state. The acceleration (deceleration) of the Universe reads 

H (t) = ~ a -Hl{t), (15) 
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implying that the type of the expansion of the Universe is exclusively determined by 7. For 7 > 2/3, the Universe 
decelerates. As shown below, in a quark-gluon plasma 7 could be as large as two times this threshold. For the viscous 
QGP model, the expressions for H will be obtained by means of a perturbative analysis and by considering numerical 
methods. 

A. Perturbative Bulk Viscosity 

In the following we consider the effects of the inclusion of the bulk viscosity in the cosmological standard model 0] 
by means of a perturbative approach. It is essential to highlight that the validity of this treatment is restricted to the 
QCD era. It is considered that the QCD era lasted a much shorter time interval than the era corresponding to a finite 
scale and a strong running coupling a s [l9l . l2(i| . As a first step in our analysis we obtain perturbatively the deviations 
of the Hubble parameter from the non-viscous cosmological picture. 

1. A First-Order Perturbative Correction 

By using a first-order perturbative correction, f(t), the solution of Eq. (|12p can be written as 

H(t)=H (t)+af(t), (17) 

where Ho(t) is the Hubble parameter corresponding to a vanishing bulk viscosity. Substituting this into Eq. (TT2"]) 
leads to 

aH (t)H (t) + |[1 + (1 - r) 7 ]aff 2 (i)ij (i) + H (t)H (t) + aH (t)f(t) 

+af(t)H (t) - (1 + r)aHl{t) + H( 7 - 2)aH$(t) + ^H 3 (t) + 9 -a~fH 2 (t)f(t) = 0. (18) 
This can be written in the form 

f{t) + G 1 {t)f{t) + G {t) = 0, (19) 

where 

G (t) = -H (t)- |[1 + (l-, ) 7]i7 oW ij„(t)-^ + (l + r)|^-|( 7 -2)i7 3 (t)-| 7 ^M, (20) 
i a ti[)(t} 4 la 

To eliminate H (t) and its derivatives, we can use Eqs. (TI~4")) . (Ti7i]) and 

H (t) = - - 3#o(t)- + 2# 3 W = (22) 
a a 07 1 6 

respectively, thus closing the coupled set of equations needed to determine the non-homogeneous coefficients Go(t) 
and G x {t). Then 

G (t) - At- 3 , (23) 
Gi{t) = 2t~\ (24) 



where the coefficient A is given by 



2 4 2 2 7-2 

A =3f [1 + (1 - r)T] ^^ (1 + r) -3V- (25) 



Equation (|19[) is a first-order linear differential equation and has the general solution given by 

f(t) = e- F ( [ e F G (t)dt + c) , (26) 



where F = j G\{t)dt, and C is an arbitrary integration constant. Substituting Eqs. (|23|) and ([24]) into Eq. ([26]) leads 
to, 

f(t) = r 2 [A\n(t)-C}. (27) 

Therefore, the solution of Eq. (fT2"|) reads 

H(t) = — - + ar 2 [Ahi(t) - C\. (28) 
37 i 

Since C = ^4 ln(io) does not depend on t, the acceleration/deceleration is obtained as 

H(t) = -JL+°[A-2(Aln(t)-C!)]. (29) 



2. A Second-Order Perturbative Correction 

A second-order perturbative correction to the Hubble parameter can be represented as 

H(t) = H (t) + af(t) + a 2 h(t). (30) 

In order to solve the resulting evolution equation for the new perturbation function h(t), we can follow the same 
procedure used in the previous Section. Substituting Eq. (|30)) and its derivatives into Eq. ([12]) reduces the problem to 
a second-order non-homogeneous differential equation given by 

f2(t)h(t) + fi(t)Ht) + fo(t)h(t) = k(t), (31) 

where the coefficients of the non-homogeneous equations are 

Ut) = \?A*&, (32) 
2 2a 3 

hit) = j(l + r)/ 2 (t)+3^[l + 7(l-r)] 1 (33) 



Aa 2 



M f ) - 3< 3 7 3 

k(t) = -a f 2 {t) 



[ 7 2 t + a (-4 + 7 (2 + 7))] (h{t)-j{l+r)h{tYj [2i-3 7 A+6 7 C + 6Aln(i)], (34) 



[-16 + 7 (8 + 7 (8 - 8r + 2,-yA + 6 7 C))] 



a 2 7 3 t 4 

3/2 (*) 
a 7 2 i 4 



[ 7 2 A(4r - 1) + 4C( 7 (3 + 7 - 2r 7 ) - 6)] 



3Aln(t)/ 2 (t) r 2 

[7 * + 2a(-6 + 7(3 + 7 - 2r 7 ))J - 

[8C*( 7 - r - 1) + 2 7 2 A(2r - 3) + 6 7 3 C*(A + C)] - (35) 



37 3 t 5 

-8 + 2 7 (2 + 7 (2-2r + 3 7 A))] - 



3 



a 2 
3 7 3 t : 

o 4 

— [ 7 2 ^(4rC* - C - rA - A) + 2C 2 ( 7 (3 + 7 - 2r 7 ) - 6)] + 
71 

4(1 + 7 -ry) 2(1 + 7 -r 7 ) 3(1 + 7 - r 7 )C 



3a s iH 3 A]n(t) Jzy ' 3a 2 7 2 t 4 Aln(t) JZW a 2 7 iMln(t) 
1 + 7 — r 7 



a 2r y 2 t 5 



h{t) [8t + &y(-A + C(4 + a)) + 3"fA{4 + a) ln(t)] 



We assume that fit), given by Eq. (f2"T|) . represents a particular solution of the following second-order homogeneous 
differential equation: 

hit) git) + hit) git) + hit) git) = 0. (36) 
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The solution g{t) of Eq. (|36| reads 



m = m 



exp(-F) 
Pit) 



dt. 



(37) 



where F = J* q f 2 (t) / ' h{t)dt. Then 



9it) 



Aln(i) - C f* t^-^l^ i\nit)Y h/fSl 



t 2 



L41n(i) - Cp 



dt. 



(38) 
(39) 



where we have denoted 



Pi = ^a 4 A, 



/3 2 = -a 4 A(l + r), 



3 H 



ry). 



A second particular solution for the homogeneous differential equation is given by 

Aln(t) - C 



t 2 



{-(2r - 3)- (1+B) r [1 + B, (2r - 3) ln(t)]} 



where T is the incomplete gamma function, and 



B = 2 



-(7-l)-l 



(40) 



(41) 



The numerical values of this second-order perturbative correction are negligibly small. Fig. [T] shows the real part of 
Eqs. ((40)) . Obviously, it is valid up to t = 1. 

The general solution of the non-homogeneous differential equation, Eq. (|3"Tj) . can be obtained with the use of the 
Wronskian determinant, W = cxp(— F), 



hit) = c 1 fit)+c 2 git)+git) / /(f) 



kit) dt 
'W)W 



fit) / 9it) 



kit) dt 
'7^)W> 



(42) 



where c\ and c 2 are arbitrary constants of integration. 

Hence, wc conclude that the first order perturbative correction hit) seems to give the dominant term, at least 
according to the quantitative comparison with the numerical solution given in the next Section. The function g(t) 
likely makes a very little contribution to the final results. Taking the non-viscous case as a particular solution leads 
to a similar conclusion. Therefore, there is no need to include the second order correction in Fig. [4] 




o.ooot 

0.70 0.75 0.80 0.85 0.90 0.95 1.00 



Fig. 1: Left panel: second particular solution git) as a function of time t. Only its real part is taken into account, 
which obviously exists only in t € [to, 1] GcV -1 . Right panel shows a comparison between first (dashed curve) and 
second (solid curve) particular solutions. 
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B. Numerical solution for the causal cosmological quark-gluon plasma 

The inclusion of the bulk viscous effects can generally be done through an effective pressure II, which is formally 
included in the effective thermodynamic pressure p e ff, where p e ff = P + n, and p is the thermodynamic pressure. 
EoS's for p and T are necessary in order to close the system of equations Eq. (fj|, ([7]) and ©. r and £ can be 
determined by using some phenomenological approaches. For instance, £ in QGP matter at high T can be estimated 
from recent lattice QCD simulations [l5l IT(|. The relaxation time can also be derived from the quasi-particle (QP) 
model, which is effectively used to reproduce the lattice QCD thermodynamics [l8[. 



P = up, 
T = (3p r , 



(43) 
(44) 



where lj — (7 — 1), u>o ~ 0.5 — 1.5 GeV, r = 0.39 and /3 = 0.718. According to recent lattice QCD simulations [15j, it 
is found that s = 1. In the quasi-particle model, it has been found that the bulk viscosity £ and the relaxation time 
r have the following barotropic forms, 



£ = (0.959 ± 0.006) p(°- 863±0 001 ), 
r = (64.189 ± 1.667) p (-o.o75±o.oos) 



(26.276 ± 1.774), 



(45) 
(46) 



which are graphically shown in Fig. [2] In the left panel, the dependence of the bulk viscosity (in units of GeV/fm 2 ) 
on p (in units of GeV/fm 3 ) is presented. The right panel shows the decay of the relaxation time r in QGP matter 
with increasing p. The parameters of the quasi-particle model are first adjusted to reproduce the thermodynamics of 
lattice QCD simulations, and then they are used to calculate both £ and r. The thermodynamic consistency of the 
barotropic relations of £ and r is guaranteed [23| . This can be partly seen from the fact that the ratio £/r strongly 
depends on p. 

Substituting Eq. (JSJ) into Eq. (J7J) leads to Eq. (TT2"]) , When assuming that the last two terms in the right hand side of 
this equation are vanishing, then we are left with a Bernoulli type differential equation, with the solutions H = and 
H fa — 27/[3a(2 — 7)]. The case H — corresponds to a vacuum Universe, and hence this solution can be ignored. If 
it would be possible to express £ as a function of t, we would obtain 2H + 3~/H 2 - 3£(<)H = 0, which is a Bernoulli 
type equation, and has the solution 



H(t) 



exp[(3/2)/g(*)rft] 
C + exp [(3/2)7 /£(*)<£]' 



(47) 



where C is an arbitrary constant of integration. A numerical solution of Eq. (TT2")) for H is presented in Fig. [5] It 
strongly depends on the initial conditions for H, H and C . As mentioned previously, the initial time is defined as 
t = 0.734 GeV" 1 , at which T = 1.0 GeV. 
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Fig. 2: The dependence of bulk viscosity £ on energy density p in QGP matter for the quasi-particle model. The right 
panel shows the p-dependence of the relaxation time r in GeV -1 from the same model. 
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IV. DISCUSSION AND FINAL REMARKS 



According to the standard cosmological standard the QCD era is characterized by a temperature range of- 0.2 < 
T < 10 GeV, or a time-interval of 18.35 < t < 0.0073 GeV _1 . In Fig. H the boundaries of the QCD era are drawn by 
vertical lines. The dashed one refers to T = 0.5 GeV, or t = 2.937 GeV -1 , respectively. The right and the left vertical 
lines give the boundaries, T = 0.2 GeV or t = 18.35 GeV" 1 and T = 1.0 GeV or t = 0.734 GeV -1 , respectively. 




Fig. 3: Various measurements of strong running coupling a s (Q) as a function of the energy scale Q in GeV units. The 
graph is taken from Ref. [l9[. The curves represent the QCD predictions, and their systematic certainty. 



In principle, the QCD era is defined according to the energy scale available in the early Universe, and according 
to the asymptotic behavior of the partonic matter. Although QCD does not predict the absolute size of the strong 
running coupling a s , it makes a precise prediction for its energy-dependence. Since a s is the most important degree 
of freedom of QCD, we define the region for which our treatment is consistent as the region where a s remains finite. 
It is important to mention that this degree of freedom appears even if the quarks are entirely excluded. The scale- 
dependence of the coupling constant a s is governed by the function/?, which encodes the running of the coupling, and 
whose perturbative expansion in the four-loop approximation is |24j . 

(3(a s (Q 2 )) = -A, a 2 s {Q 2 ) - ft a 3 s (Q 2 ) - ft c^(Q 2 ) - ft a 5 s (Q 2 ) + • • • , (48) 

where the coefficients ft, ft, ft, and ft have been determined by perturbative methods, and found to be dependent 
on the number of active quark flavors n/ at the corresponding energy scale Q. For simplicity, let us consider the 
one-loop approximation. Then 

asiQ2) ~ "(33 - 271/) ln(Q 2 /A 2 )' (49) 

where A is the QCD energy scale, depending onn/, and on the renormalization scheme. If Q = A, then a s diverges. 
On the other hand, a s — > at very large energy scale. A summary of a s -measurements is graphically displayed in 
Fig- El It is obvious that the measurements unambiguously confirm the QCD predictions for a s , i.e, Eq. (|49|) . as well 
as the approach towards asymptotic freedom. The figure is taken from Ref. The QCD era is conjectured to 

start from a very high energy scale, i.e, very small a s , and ends up when Q — > A, i.e, at divergent a s . As a matter 
of precaution, we assume that the QGP matter exists in a much narrower energy scale, up to just 3 — 5 times A. 
According to the recent lattice QCD simulations fl5| . A can be localized at T c ~ 0.2 GeV, which sets the latest end of 
the QCD era. Also, we know that the bulk viscosity is about to diverge at T Cl and decays with increasing temperature. 
Based on these two ingredients, we set some narrow limits to the validity of including finite bulk viscosity in the QGP 
era. 

In the left panel of Fig. 01 the evolution of the Hubble parameter H is given for different values of the first-order 
bulk viscosity coefficient a. Obviously, the deviation from the non-viscous evolution increases with increasing a. The 
vertical lines determine when or where QGP matter is becoming dominant in the early Universe. Pre- and post- 
cosmological eras are probably characterized by EoS's that might be different from the ones used in the present work. 
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t [1/GeV] t [1/GeV] 

Fig. 4: Left panel: the time evolution of the Hubble parameter H(t) is depicted for different values of the first-order 
perturbative bulk viscosity, a. The vertical lines determine the time t and temperature T boundaries of the QCD era 
in early Universe that have been suggested in present work. In the right panel, H (t) is given for different values of a. 

H (time) 




Fig. 5: Numerical solution for the Hubble parameter H(t) for a bulk viscous fluid with finite bulk viscosity in the 
full causal Israel- Stewart theory (dashed-dotted bottom curve). The dashed curve (top) represents H(t) in the QCD 
era with non-viscous background matter. The physical units in both axis are GeV. The solid curve (left) gives the 
evolution in the pre-QCD eras, where vanishing viscosity characterizes the background matter. 

Under these assumptions, the viscosity effects are assumed to be localized within the QGP era. In the right panel, 
the evolution of the time derivative of H is presented for different values of a. As shown in Eq. (TT4")) , H is the sum 
of two terms. The first term gives a ratio of acceleration, a, to a scale parameter, a. The second term is the square 
of H itself. The sum of these two terms is negative at small t. At large t, it the two terms are identical not only 
qualitatively, but also quantitatively That the Universe is obviously accelerating, a > 0, is also shown by the values 
of the scale parameter a, which increase with increasing t. Therefore, we conclude that including finite bulk viscosity 
seems to affect the expansion of the Universe in a significant way. The evolution of the Hubble parameter likely decays 
with a fast rate, that increases with increasing t. Furthermore, we point out that increasing the values of the bulk 
viscosity coefficient a accelerates the decay of H . 

A numerical solution for the evolution equation, Eq. (TT2)) , at finite bulk viscosity, is presented in Fig. [SJ Here, we 
take into consideration the barotropic dependence of £, as it has been obtained in the lattice QCD simulations. We do 
not make any further assumptions on it. Also, the dependence of the relaxation time r on the energy density has not 
been modified. These two ingredients have been strictly implemented. The consistency of these barotropic EoS's with 
the laws of thermodynamics has been discussed in Ref. [23|. It seems that all barotropic EoS's used in the present (and 
also in the previous work 0, [HJ) are thermodynamically consistent. Therefore, we numerically solve the evolution 
equation for the Hubble parameter. The analytic solution has been obtained in Ref. 0- This solution is sensitive 
to the particular assumptions made in order to solve the second Abel-type non-linear non-homogeneous differential 
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equation. To avoid their effects, we consider the numerical solution to the evolution equation. The boundary conditions 
required for the numerical methods are defined at the boundaries of QGP era in the early Universe. 

A comparison between the time evolution of the Hubble parameter H in the non-viscous and viscous background 
matter is presented in Fig. [5] The top dashed curve represents the evolution of the non-viscous cosmological matter, 
given by Eq. (TH|) . The dotted-dashed curve gives the evolution of the viscous cosmological matter. In order to 
compare with the results presented in Fig. |4] the earliest (left) time boundary has been zoomed in. By comparing the 
numerical and analytical results leads us to the conclusion that the two methods (perturbative and numerical) agree 
in reproducing the evolution of H. Also, we conclude that the time evolution of H for viscous cosmological matter 
is faster than the evolution for non-viscous cosmological matter. Such a difference can be estimated, quantitatively 
Therefore, essential cosmological consequences are to be expected from the inclusion of bulk viscous effects in the 
description of the cosmological expansion of the Universe. 
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